#include "Mesh.h"

#include <cmath>
#include <cfloat>
#include <cassert>
#include <algorithm>
#include <fstream>

#define ACCURACY 1e-7

using namespace Mesh;

// 基本运算
vector3f vector3f::operator*(double t) const
{
	return {x * t, y * t, z * t};
}

vector3f vector3f::operator/(double t) const
{
	assert(fabs(t) > ACCURACY);
	return {x / t, y / t, z / t};
}

vector3f vector3f::operator-(const vector3f &v) const
{
	return {x - v.x, y - v.y, z - v.z};
}

vector3f vector3f::operator+(const vector3f &v) const
{
	return {x + v.x, y + v.y, z + v.z};
}

vector3f Mesh::operator-(const vector3f &v)
{
	return {-v.x, -v.y, -v.z};
}

vector3f point3f::operator-(const point3f &p) const
{
	return {x - p.x, y - p.y, z - p.z};
}

// 顶点 + 向量 = 顶点
point3f Mesh::operator+(const point3f &p, const vector3f &v)
{
	return {p.x + v.x, p.y + v.y, p.z + v.z};
}

// 向量 + 顶点 = 顶点
point3f Mesh::operator+(const vector3f &v, const point3f &p)
{
	return {p.x + v.x, p.y + v.y, p.z + v.z};
}

// 获得绕轴 v 旋转 theta 角的四元数
quaternion::quaternion(vector3f v, double theta) : s(cos(theta / 2)), v(normalize(v) * sin(theta / 2)) {}

// 四元数基本运算
quaternion quaternion::operator*(double t) const
{
	return quaternion(s * t, v * t);
}
quaternion quaternion::operator/(double t) const
{
	assert(fabs(t) > ACCURACY);
	return quaternion(s / t, v / t);
}

quaternion quaternion::operator+(const quaternion &q) const
{
	return quaternion(s + q.s, v + q.v);
}
quaternion quaternion::operator-(const quaternion &q) const
{
	return quaternion(s - q.s, v - q.v);
}

quaternion quaternion::operator*(const quaternion &q) const
{
	return quaternion(s * q.s - dot(v, q.v), q.v * s + v * q.s + cross(v, q.v));
}
quaternion quaternion::operator/(const quaternion &q) const
{
	// 先获得 q 的共轭 / |q|^2，这就是 q^-1，然后直接调用乘法
	return (*this) * quaternion(q.s, -q.v) / pow(length(q), 2);
}

// 基本运算
coord4f coord4f::operator*(double t) const
{
	coord4f c = *this;
	if (w == 1)
	{
		c.x /= t;
		c.y /= t;
		c.z /= t;
	}
	else
	{
		c.x *= t;
		c.y *= t;
		c.z *= t;
	}
	return c;
}

coord4f coord4f::operator/(double t) const
{
	assert(fabs(t) > ACCURACY);
	coord4f c = *this;
	if (w == 1)
	{
		c.x *= t;
		c.y *= t;
		c.z *= t;
	}
	else
	{
		c.x /= t;
		c.y /= t;
		c.z /= t;
	}
	return c;
}

coord4f coord4f::operator-(const coord4f &v) const
{
	return {x - v.x, y - v.y, z - v.z, w - v.w};
}

coord4f coord4f::operator+(const coord4f &v) const
{
	coord4f c = {x + v.x, y + v.y, z + v.z, w + v.w};
	if (c.w == 2)
	{
		c.x /= c.w;
		c.y /= c.w;
		c.z /= c.w;
	}
	return c;
}

coord4f Mesh::operator-(const coord4f &v)
{
	return {-v.x, -v.y, -v.z, v.w};
}

// 颜色直接叠加
rgb_space rgb_space::operator+(const rgb_space &c1) const
{
	rgb_space c = {c1.r + r, c1.g + g, c1.b + b};
	c.r = (c.r > 1) ? 1 : c.r;
	c.g = (c.g > 1) ? 1 : c.g;
	c.b = (c.b > 1) ? 1 : c.b;
	return c;
}

// 最后释放内存
HEMesh::~HEMesh()
{
	for (auto it = m_vList.begin(); it != m_vList.end(); ++it)
		if ((*it) != nullptr)
			delete *it;
	for (auto it = m_eList.begin(); it != m_eList.end(); ++it)
		if ((*it) != nullptr)
			delete *it;
	for (auto it = m_fList.begin(); it != m_fList.end(); ++it)
		if ((*it) != nullptr)
			delete *it;
}

// 创建顶点，并返回指针
vertex3f *HEMesh::CreateVertex(point3f p)
{
	vertex3f *vp = new vertex3f({p, nullptr});
	m_vList.push_back(vp);
	return vp;
}

// 创建三角面，并自动衔接
face *HEMesh::CreateFace(vertex3f *v0, vertex3f *v1, vertex3f *v2)
{
	// 判断三点是否共线
	if (onLine(v0->pos, v1->pos, v2->pos))
		return nullptr;

	face *f = new face;

	// 创建三条半边
	half_edge *e = new half_edge;
	half_edge *ne = new half_edge;
	half_edge *pe = new half_edge;
	m_eList.push_back(e);
	m_eList.push_back(ne);
	m_eList.push_back(pe);

	// 面持有任意一个半边
	f->e = e;

	// 初始化半边
	e->v = v0;
	e->incFace = f;
	e->next = ne;
	e->prev = pe;
	e->twin = nullptr;

	ne->v = v1;
	ne->incFace = f;
	ne->next = pe;
	ne->prev = e;
	ne->twin = nullptr;

	pe->v = v2;
	pe->incFace = f;
	pe->next = e;
	pe->prev = ne;
	pe->twin = nullptr;

	// 遍历整个面表，逐一尝试建立 twin 连接
	// 优化：不再遍历面表，而是选择与顶点相邻的面尝试建立 twin 连接；当面的数量较大时具有优势
	// 再优化：由于三角形顶点两两共用面，只需要计算两个顶点的邻面
	// 注意要先取出所有邻面再建立连接，防止出现与自己连接的情况
	std::deque<face *> fList0 = adjacentFace(v0);
	std::deque<face *> fList1 = adjacentFace(v1);
	for (auto it = fList0.begin(); it != fList0.end(); ++it)
		stitchFace(f, *it);
	for (auto it = fList1.begin(); it != fList1.end(); ++it)
		stitchFace(f, *it);

	// 如果顶点没有半边，就赋予半边，放在这里，因为如果顶点没有半边，不需要计算邻面
	if (v0->e == nullptr)
		v0->e = e;
	if (v1->e == nullptr)
		v1->e = ne;
	if (v2->e == nullptr)
		v2->e = pe;

	// 最后推入面
	m_fList.push_back(f);

	return f;
}

// 创建两个三角面并返回，自动判定衔接相邻面
face *HEMesh::CreateFace(vertex3f *v0, vertex3f *v1, vertex3f *v2, vertex3f *v3)
{
	CreateFace(v0, v1, v3);
	CreateFace(v1, v2, v3);
}

// 删除顶点
void HEMesh::DeleteVertex(vertex3f *v)
{
	// 获得所有相邻的面，将它们删除
	std::deque<face *> fList = adjacentFace(v);
	for (auto it = fList.begin(); it != fList.end(); ++it)
		DeleteFace(*it);

	// 最后删除顶点
	m_vList.erase(std::find(m_vList.begin(), m_vList.end(), v));
	delete v;
	v = nullptr;
}

// 删除面
void HEMesh::DeleteFace(face *f)
{
	// 获得三条半边
	half_edge *e[] = {f->e, f->e->next, f->e->prev};

	// 调整每条半边起点保存的半边
	for (int i = 0; i < 3; i++)
	{
		// 获得半边的起源点，和以 v 为起点的半边
		vertex3f *v = e[i]->v;

		// 当起点保存的半边是此半边时，就修改保存的半边。因此可以提前修改，这里就不会修改
		if (v->e == e[i])
		{
			// 删除面时，如果面上的顶点已经孤立，说明出现了非流形网格
			// 此时邻边的数量为零，应当重点注意
			std::deque<half_edge *> eList = adjacentHalfEdge(v);

			// 如果 v 不止一条半边，就赋予新的半边，否则 v 将会孤立
			for (auto it = eList.begin(); it != eList.end(); ++it)
			{
				// 赋予不同的半边
				if ((*it) != e[i])
				{
					v->e = *it;
					break;
				}
			}
		}
	}

	// 删除每条半边
	for (int i = 0; i < 3; i++)
	{
		// 首先将每条半边的反向边的反向边设为空
		half_edge *te = e[i]->twin;
		if (te != nullptr)
			te->twin = nullptr;

		// 如果此半边的起点仍然持有此半边，那么起点将会孤立，将该点保有的半边设为空
		if (e[i]->v->e == e[i])
		{
			e[i]->v->e = nullptr;
			e[i]->v = nullptr;
		}

		// 然后就可以删除此半边
		m_eList.erase(std::find(m_eList.begin(), m_eList.end(), e[i]));
		delete e[i];
		e[i] = nullptr;
	}

	// 最后删除面
	m_fList.erase(std::find(m_fList.begin(), m_fList.end(), f));
	delete f;
	f = nullptr;
}

// 清除孤立点
void HEMesh::DeleteDiscrete()
{
	// 找到所有孤立点
	for (auto it = m_vList.begin(); it != m_vList.end();)
	{
		// 如果不存在以它为起点的半边，说明此顶点孤立，直接删除
		if ((*it)->e == nullptr)
			it = m_vList.erase(it);
		else
			++it;
	}
}

// 是否有法向
bool HEMesh::HasVertexNormal() const
{
	return !m_nList.empty();
}

// 是否有颜色
bool HEMesh::HasVertexColor() const
{
	return !m_cList.empty();
}

// 判断网格是否是凸包
bool HEMesh::IsConvexHull() const
{
	// 对每个顶点，遍历每一个面
	for (auto vit = m_vList.begin(); vit != m_vList.end(); ++vit)
	{
		// 要排除离散点
		if ((*vit)->e == nullptr)
			continue;

		// 判断顶点是不是在所有面的上方，如果在某个面的上面，说明不是凸包
		for (auto fit = m_fList.begin(); fit != m_fList.end(); ++fit)
			if (aboveFace((*vit)->pos, *fit))
				return false;
	}

	return true;
}

// 判断点是否在网格的凸包中
bool HEMesh::IsInMesh(point3f p) const
{
	// 遍历每一个面，p 如果在某个面的上面，说明 p 在网格外
	for (auto it = m_fList.begin(); it != m_fList.end(); ++it)
		if (aboveFace(p, *it))
			return false;
	return true;
}

// 颜色线性插值
rgb_space Mesh::linear(rgb_space c1, rgb_space c2, double t)
{
	return {c1.r * (1 - t) + c2.r * t,
			c1.g * (1 - t) + c2.g * t,
			c1.b * (1 - t) + c2.b * t};
}

// 点线性插值
point3f Mesh::linear(point3f p1, point3f p2, double t, double a, double b)
{
	return {(p1.x * (b - t) + p2.x * (t - a)) / (b - a),
			(p1.y * (b - t) + p2.y * (t - a)) / (b - a),
			(p1.z * (b - t) + p2.z * (t - a)) / (b - a)};
}

// 计算点的距离
double Mesh::distance(point3f p1, point3f p2)
{
	return sqrt(pow(p1.x - p2.x, 2) + pow(p1.y - p2.y, 2) + pow(p1.z - p2.z, 2));
}

// 行列式 (v1 x v2) * v3
double Mesh::det(vector3f v1, vector3f v2, vector3f v3)
{
	return dot(cross(v1, v2), v3);
}

// 向量点积
double Mesh::dot(vector3f v1, vector3f v2)
{
	return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
}

// 向量叉积
vector3f Mesh::cross(vector3f v1, vector3f v2)
{
	return {v1.y * v2.z - v1.z * v2.y,
			-v1.x * v2.z + v1.z * v2.x,
			v1.x * v2.y - v1.y * v2.x};
}

// 向量规范化
vector3f Mesh::normalize(vector3f p)
{
	return p / length(p);
}
// 四元数规范化
quaternion Mesh::normalize(quaternion q)
{
	return q / length(q);
}

// 向量长度
double Mesh::length(vector3f p)
{
	return sqrt(p.x * p.x + p.y * p.y + p.z * p.z);
}
// 四元数长度
double Mesh::length(quaternion q)
{
	return sqrt(q.s * q.s + pow(length(q.v), 2));
}

// 判断三点是否共线
bool Mesh::onLine(point3f p1, point3f p2, point3f p3)
{
	vector3f vec1 = p1 - p2;
	vector3f vec2 = p1 - p3;

	// 叉积足够大，则认为不共线
	if (length(cross(vec1, vec2)) > ACCURACY)
		return false;

	return true;
}

// 获得从轴 a 转动到轴 b 的四元数
quaternion Mesh::rotate(vector3f a, vector3f b)
{
	// 确保两个轴长度不为零
	assert(length(a) > ACCURACY && length(b) > ACCURACY);

	// 先把 a 转动到 z 轴：先绕 z 轴转，然后绕 y 轴转
	quaternion q1(1), q2(1);
	double proj = sqrt(pow(a.x, 2) + pow(a.y, 2));
	// 如果已经在 z 轴，就不动；否则计算两个四元数
	if (proj > ACCURACY)
	{
		q1 = quaternion({0, 0, 1}, -acos(a.x / proj));
		q2 = quaternion({0, 1, 0}, -acos(a.z / length(a)));
	}

	// 然后把 b 同样转动到 z 轴
	quaternion q3(1), q4(1);
	proj = sqrt(pow(b.x, 2) + pow(b.y, 2));
	// 如果已经在 z 轴，就不动；否则计算两个四元数
	if (proj > ACCURACY)
	{
		q3 = quaternion({0, 0, 1}, -acos(b.x / proj));
		q4 = quaternion({0, 1, 0}, -acos(b.z / length(b)));
	}
	return quaternion(1) / q3 / q4 * q2 * q1;
}

// “缝合”两个面，返回是否成功
bool Mesh::stitchFace(face *f1, face *f2)
{
	// 先取得公共半边
	std::deque<half_edge *> e = commonOnFace(f1, f2);

	// 如果有公共边就连接
	if (e[0] != nullptr && e[1] != nullptr)
	{
		// 如果公共边已经建立，但是不是两个面的公共边，则返回 false
		if ((e[0]->twin == nullptr && e[1]->twin == nullptr) ||
			(e[0]->twin->twin == e[0] && e[1]->twin->twin == e[1]))
		{
			// 然后互相为 twin
			e[0]->twin = e[1];
			e[1]->twin = e[0];
			return true;
		}
	}
	// 报错说明这个面破坏了网格的流形结构，或者没有公共边
	return false;
}

// 判断一个点是否在面上
bool Mesh::onFace(point3f p, face *f)
{
	// 排除面是空的情况
	if (f == nullptr)
		return false;

	// 找到面的一个顶点
	point3f p0 = f->e->v->pos;
	vector3f v = p - p0;

	// 如果与法向正交，说明共面
	return (fabs(dot(v, faceNormal(f))) < ACCURACY) ? true : false;
}

// 判断一个点是否在面的上方
bool Mesh::aboveFace(point3f p, face *f)
{
	// 排除面是空的情况
	if (f == nullptr)
		return false;

	// 取一个顶点计算夹角
	point3f p0 = f->e->v->pos;
	vector3f v = p - p0;

	// 这里精度需要仔细考虑，不然判断会有问题
	return (dot(v, faceNormal(f)) > ACCURACY) ? true : false;
}

// Gouraud 插值：计算顶点的平均法向
vector3f Mesh::Gouraud(vertex3f *v)
{
	// 先获取所有相邻的面
	std::deque<face *> fList = adjacentFace(v);

	// 累加然后规范化
	vector3f n = {0, 0, 0};
	for (auto it = fList.begin(); it != fList.end(); ++it)
		n = n + faceNormal(*it);
	return normalize(n);
}

// 计算单位面法向
vector3f Mesh::faceNormal(face *f)
{
	// 找到面的三个顶点计算法向
	point3f v0 = f->e->v->pos;
	point3f v1 = f->e->next->v->pos;
	point3f v2 = f->e->next->next->v->pos;
	return normalize(cross(v1 - v0, v2 - v1));
}

// 获得边的目标点
vertex3f *Mesh::destination(half_edge *e)
{
	if (e == nullptr)
		return nullptr;
	e = e->next;
	return e->v;
}

// 获得一个面上的所有顶点
std::deque<vertex3f *> Mesh::vertexOnFace(face *f)
{
	std::deque<vertex3f *> vList;
	half_edge *e = f->e;
	for (int i = 0; i < 3; i++)
	{
		vList.push_back(e->v);
		e = e->next;
	}
	return vList;
}

// 获得一个面上的所有半边
std::deque<half_edge *> Mesh::halfEdgeOnFace(face *f)
{
	std::deque<half_edge *> eList;
	half_edge *e = f->e;
	for (int i = 0; i < 3; i++)
	{
		eList.push_back(e);
		e = e->next;
	}
	return eList;
}

// 获取两个面的公共半边
std::deque<half_edge *> Mesh::commonOnFace(face *f0, face *f1)
{
	std::deque<half_edge *> eList(2, nullptr);

	half_edge *e0 = f0->e;
	half_edge *e1 = f1->e;

	// 遍历两个面上的所有边
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			// e0 的起点是 e1 的终点，e1 的起点是 e0 的终点时，说明是公共边，储存后直接返回
			if (e0->v == destination(e1) && e1->v == destination(e0))
			{
				eList[0] = e0;
				eList[1] = e1;
				return eList;
			}
			e1 = e1->next;
		}
		// 注意 e1 返回原位
		e1 = f1->e;
		e0 = e0->next;
	}
	return eList;
}

// 获得相邻的所有顶点
std::deque<vertex3f *> Mesh::adjacentVertex(vertex3f *v)
{
	std::deque<vertex3f *> vList;

	// 先获得以 v 为起点的半边，其目标点就是相邻顶点
	std::deque<half_edge *> eList = adjacentHalfEdge(v);
	for (int i = 0; i < eList.size(); i++)
		vList.push_back(destination(eList[i]));

	return vList;
}

// 获得以 v 为起点的半边
std::deque<half_edge *> Mesh::adjacentHalfEdge(vertex3f *v)
{
	std::deque<half_edge *> eList;
	half_edge *e = v->e;

	// 先顺时针寻找
	while (e != nullptr)
	{
		eList.push_back(e);
		e = e->twin;

		// 反向边，然后找到下一条边
		if (e != nullptr)
			e = e->next;

		// 如果回到初始边，就退出
		if (e == v->e)
			break;
	}

	// 如果没有回到初始边，就逆时针寻找，排除孤立顶点
	if (e == nullptr && v->e != nullptr)
	{
		e = v->e->prev->twin;
		while (e != nullptr)
		{
			eList.push_back(e);
			e = e->prev->twin;
		}
		// 这里不需要判断回到初始边，因为如果一个方向不能返回，那么反过来也不能返回
	}
	return eList;
}

// 获得相邻的所有面
std::deque<face *> Mesh::adjacentFace(vertex3f *v)
{
	// 先获得以 v 为起点的半边
	std::deque<half_edge *> eList = adjacentHalfEdge(v);

	// 每条半边左边的面就是相邻的面
	std::deque<face *> fList;
	for (auto it = eList.begin(); it != eList.end(); ++it)
		fList.push_back((*it)->incFace);

	return fList;
}

// 获得相邻的所有面
std::deque<face *> Mesh::adjacentFace(face *f)
{
	std::deque<face *> fList;
	half_edge *e = f->e;
	for (int i = 0; i < 3; i++)
	{
		if (e->twin != nullptr)
			fList.push_back(e->twin->incFace);
		e = e->next;
	}
	return fList;
}

// 尝试用网格点创建四面体，返回是否成功
bool Mesh::tetrahedron(HEMesh *mesh)
{
	// 总顶点数至少是 4，并且还没有创建任何面和边
	if (mesh == nullptr || mesh->m_vList.size() < 4 || mesh->m_fList.size() > 0 || mesh->m_eList.size() > 0)
		return false;

	// 先用前 3 个顶点建立一个三角形
	auto it = mesh->m_vList.begin();
	vertex3f *v1 = *it++;
	vertex3f *v2 = *it++;
	vertex3f *v3 = nullptr;
	vertex3f *v4 = nullptr;

	// 第三个顶点不能与前两个共线
	for (; it != mesh->m_vList.end(); ++it)
	{
		// 找到不共线的顶点，方便起见，将它插入到头部
		if (!onLine(v1->pos, v2->pos, (*it)->pos))
		{
			v3 = *it;
			mesh->m_vList.erase(it);
			mesh->m_vList.push_front(v3);
			break;
		}
	}
	// 如果没找到，说明全共线，退出
	if (v3 == nullptr)
		return false;

	// 创建三角形
	face *f = mesh->CreateFace(v1, v2, v3);

	// 寻找第 4 个顶点，不与前三个共面
	it = mesh->m_vList.begin();
	it++, it++, it++;
	for (; it != mesh->m_vList.end(); ++it)
	{
		// 找到不共面的顶点，方便起见，将它插入到头部
		if (!onFace((*it)->pos, f))
		{
			v4 = *it;
			mesh->m_vList.erase(it);
			mesh->m_vList.push_front(v4);
			break;
		}
	}
	// 如果没找到，说明全共面，退出
	if (v4 == nullptr)
		return false;

	// 检查 4 顶点是否在面的上方
	bool above = aboveFace(v4->pos, f);
	if (above)
	{
		// 如果在上方，就要删除原来的面，将顺序反转，以确保法向正确
		mesh->DeleteFace(f);
		mesh->CreateFace(v3, v2, v1);
		mesh->CreateFace(v1, v4, v3);
		mesh->CreateFace(v3, v4, v2);
		mesh->CreateFace(v1, v2, v4);
	}
	else
	{
		mesh->CreateFace(v1, v3, v4);
		mesh->CreateFace(v3, v2, v4);
		mesh->CreateFace(v1, v4, v2);
	}

	return true;
}

// 取网格点的凸包
bool Mesh::convexHull(HEMesh *mesh)
{
	// 先尝试创建一个四面体
	if (!tetrahedron(mesh))
		return false;

	// 增量算法，对每个四面体以外的顶点，计算被这个顶点照射到的面
	// 将这些面删除，然后把面的边与顶点建立新的面
	auto vit = mesh->m_vList.begin();
	vit++, vit++, vit++, vit++;
	for (; vit != mesh->m_vList.end(); ++vit)
	{
		// 标记被照射的面，只需要保存一个面的指针，因为所有标记面都相邻
		std::deque<face *> fList;
		for (auto fit = mesh->m_fList.begin(); fit != mesh->m_fList.end(); ++fit)
		{
			// 在表面上面，说明被照射到，将其记录下来
			if (aboveFace((*vit)->pos, *fit))
			{
				(*fit)->identity = 1;
				fList.push_back(*fit);
			}
			else
				(*fit)->identity = 0;
		}

		// 检查被照射的面，标记相邻的未照射面
		std::deque<half_edge *> eList;
		for (auto fit = fList.begin(); fit != fList.end(); fit++)
		{
			half_edge *e = (*fit)->e;
			// 这里不需要考虑没有反向边，因为凸包的边一定有反向边
			for (int i = 0; i < 3; i++)
			{
				// 记录边界边，同时将边界边上的顶点保存的半边修改为对边，防止在删除面时孤立
				if (e->twin->incFace->identity != 1)
				{
					destination(e)->e = e->twin;
					eList.push_back(e->twin);
				}
				e = e->next;
			}
		}

		// 删除光照面。直接删除不会有问题，如果出现顶点指向的边为空的情况
		// 有两种可能：1. 离散点，则在最后会被删除；2. 处于被光照和未被光照的边界上，这种情况在前面进行修正
		for (auto it = fList.begin(); it != fList.end(); ++it)
			mesh->DeleteFace(*it);

		// 遍历所有反向边，按照连接顺序，每次要创建具有公共边的新面，防止出现非流形网格
		// 这样创建的面都会相邻，不会出现恰好有一个公共顶点的两个面，而周围都是空的情况
		int count = eList.size();
		auto eit = eList.begin();
		while (count > 0)
		{
			mesh->CreateFace(destination((*eit)), (*eit)->v, *vit);
			count--;

			// 找到目标点对应的边
			for (auto it = eList.begin(); it != eList.end(); ++it)
			{
				if ((*it)->v == destination((*eit)))
				{
					eit = it;
					break;
				}
			}
		}
	}
	// 最后移除所有离散点
	mesh->DeleteDiscrete();
	return true;
}

// 修复网格，有时删除操作造成某些非离散点存放的半边为空，此函数修复这一问题，返回是否进行了修复
bool Mesh::fixMesh(HEMesh *mesh)
{
	bool isBad = false;
	for (auto it = mesh->m_fList.begin(); it != mesh->m_fList.end(); ++it)
	{
		// 获得三条半边
		half_edge *e[] = {(*it)->e, (*it)->e->next, (*it)->e->prev};
		for (int i = 0; i < 3; i++)
		{
			if (e[i]->v->e == nullptr)
			{
				e[i]->v->e = e[i];
				isBad = true;
			}
		}
	}
	return isBad;
}

// 球面：中心位置、三个半长轴，分割份数
bool Mesh::sphere(HEMesh *mesh, point3f center, double a, double b, double c, int div)
{
	// 首先检查参数有效性，分割份数至少应该是 2
	if (mesh == nullptr || div < 1)
		return false;

	// 把 a,b,c 都换成正的
	a = fabs(a);
	b = fabs(b);
	c = fabs(c);

	// 使用极坐标取点，获得纵向步长，也是横向步长
	const double PI = 3.1415926;
	double dphi = PI / div;
	double phi = dphi;
	double theta = 0;

	// 不创建极点，i 是纵向，因此纵向点的数量少一些，j 是横向，横向分割份数是纵向的 2 倍
	for (int i = 0; i < div - 1; i++)
	{
		for (int j = 0; j < div * 2; j++)
		{
			double x = a * sin(phi) * cos(theta) + center.x;
			double y = b * sin(phi) * sin(theta) + center.y;
			double z = c * cos(phi) + center.z;
			mesh->CreateVertex({x, y, z});
			theta += dphi;
		}
		theta = 0;
		phi += dphi;
	}

	// 现在创建这些点构成的面，每一行和上一行构成面
	for (int i = 1; i < div - 1; i++)
	{
		for (int j = 0; j < div * 2; j++)
		{
			vertex3f *v0 = mesh->m_vList[i * (div * 2) + j];
			vertex3f *v1 = mesh->m_vList[i * (div * 2) + (j + 1) % (div * 2)];
			vertex3f *v2 = mesh->m_vList[(i - 1) * (div * 2) + (j + 1) % (div * 2)];
			vertex3f *v3 = mesh->m_vList[(i - 1) * (div * 2) + j];
			mesh->CreateFace(v0, v1, v2, v3);
		}
	}

	// 最后创建极点，连接极点和相邻的点
	mesh->CreateVertex({center.x, center.y, center.z + c});
	for (int j = 0; j < div * 2; j++)
	{
		vertex3f *v0 = mesh->m_vList[j];
		vertex3f *v1 = mesh->m_vList[(j + 1) % (div * 2)];
		vertex3f *v2 = mesh->m_vList.back();
		mesh->CreateFace(v0, v1, v2);
	}

	mesh->CreateVertex({center.x, center.y, center.z - c});
	for (int j = 0; j < div * 2; j++)
	{
		vertex3f *v0 = mesh->m_vList[(div - 2) * (div * 2) + j];
		vertex3f *v1 = mesh->m_vList[(div - 2) * (div * 2) + (j + 1) % (div * 2)];
		vertex3f *v2 = mesh->m_vList.back();
		mesh->CreateFace(v1, v0, v2);
	}
	return true;
}

// 柱面：中心位置，中轴，半径，高度，分割份数
bool Mesh::cylinder(HEMesh *mesh, point3f center, vector3f axis, double r, double h, int divr, int divh)
{
	// 首先检查参数有效性，分割份数至少应该是 2
	if (mesh == nullptr || divr < 1 || divh < 1)
		return false;

	// 把 r,h 都换成正的
	r = fabs(r);
	h = fabs(h);

	// 使用极坐标取点，获得纵向步长，和横向步长
	const double PI = 3.1415926;
	double dtheta = 2 * PI / divr;
	double theta = 0;
	double dh = h / divh;

	// 获得从 z 轴转动到 axis 的四元数
	quaternion q = rotate({0, 0, 1}, axis);

	double nz = -h / 2;
	// 从下向上依次创建顶点
	for (int i = 0; i < divh + 1; i++)
	{
		for (int j = 0; j < divr; j++)
		{
			// 先记录一个纯四元数，将其转动
			quaternion p({r * cos(theta), r * sin(theta), nz});
			p = q * p / q;
			mesh->CreateVertex({p.v.x + center.x, p.v.y + center.y, p.v.z + center.z});
			theta += dtheta;
		}
		theta = 0;
		nz += dh;
	}

	// 现在创建这些点构成的面，每一行和上一行构成面
	for (int i = 0; i < divh; i++)
	{
		for (int j = 0; j < divr; j++)
		{
			vertex3f *v0 = mesh->m_vList[i * divr + j];
			vertex3f *v1 = mesh->m_vList[i * divr + (j + 1) % divr];
			vertex3f *v2 = mesh->m_vList[(i + 1) * divr + (j + 1) % divr];
			vertex3f *v3 = mesh->m_vList[(i + 1) * divr + j];
			mesh->CreateFace(v0, v1, v2, v3);
		}
	}
	return true;
}

// 网格保存
bool Mesh::saveMesh(const std::string &filename, HEMesh *mesh)
{
	if (mesh == nullptr)
		return false;

	// 首先输出各个属性的描述
	std::fstream fp(filename, std::ios::out);
	fp << "VertexNum " << mesh->VertexNum() << std::endl;
	fp << "FaceNum " << mesh->FaceNum() << std::endl;
	fp << "HasColor " << mesh->HasVertexColor() << std::endl;
	fp << "HasNormal " << mesh->HasVertexNormal() << std::endl;

	// 然后输出顶点信息，并对顶点标号
	for (int i = 0; i < mesh->VertexNum(); i++)
	{
		point3f &pos = mesh->m_vList[i]->pos;
		mesh->m_vList[i]->identity = i;
		fp << "v " << pos.x << " " << pos.y << " " << pos.z << std::endl;
	}

	// 输出面信息
	for (int i = 0; i < mesh->FaceNum(); i++)
	{
		fp << "f ";
		std::deque<vertex3f *> vList = vertexOnFace(mesh->m_fList[i]);
		for (int j = 0; j < 3; j++)
			fp << vList[j]->identity << " ";
		fp << std::endl;
	}

	// 如果有颜色就输出颜色
	if (mesh->HasVertexColor())
	{
		for (int i = 0; i < mesh->VertexNum(); i++)
		{
			rgb_space &rgb = mesh->m_cList[i];
			fp << "c " << rgb.r << " " << rgb.g << " " << rgb.b << std::endl;
		}
	}

	// 如果有顶点法向就输出顶点法向
	if (mesh->HasVertexNormal())
	{
		for (int i = 0; i < mesh->VertexNum(); i++)
		{
			vector3f &n = mesh->m_nList[i];
			fp << "n " << n.x << " " << n.y << " " << n.z << std::endl;
		}
	}

	// 最后输出终止字符
	fp << "finish" << std::endl;

	return true;
}

// 网格读取
bool Mesh::loadMesh(const std::string &filename, HEMesh *mesh)
{
	if (mesh == nullptr)
		return false;

	// 文件打不开报错
	std::fstream fp(filename, std::ios::in);
	if (!fp.is_open())
		return false;

	std::string detail;
	double value;

	// 首先跳过各个属性的描述
	fp >> detail >> value;
	fp >> detail >> value;
	fp >> detail >> value;
	fp >> detail >> value;

	// 然后根据第一个字符判断读取内容
	fp >> detail;
	while (detail == "v")
	{
		double x, y, z;
		fp >> x >> y >> z;
		mesh->CreateVertex({x, y, z});
		fp >> detail;
	}
	while (detail == "f")
	{
		int i1, i2, i3;
		fp >> i1 >> i2 >> i3;
		vertex3f *v1 = mesh->m_vList[i1];
		vertex3f *v2 = mesh->m_vList[i2];
		vertex3f *v3 = mesh->m_vList[i3];
		mesh->CreateFace(v1, v2, v3);
		fp >> detail;
	}
	while (detail == "c")
	{
		double r, g, b;
		fp >> r >> g >> b;
		mesh->m_cList.push_back(rgb_space(r, g, b));
		fp >> detail;
	}
	while (detail == "n")
	{
		double x, y, z;
		fp >> x >> y >> z;
		mesh->m_nList.push_back(vector3f(x, y, z));
		fp >> detail;
	}

	return true;
}